Sound propagation in inhomogeneous waveguides with sound-speed profiles using the multimodal admittance method
Li Qi1, 2, 3, Liu Juan1, 2, 3, †, Guo Wei1, 2, 3
Acoustic Science and Technology Laboratory, Harbin Engineering University, Harbin 150001, China
Key Laboratory of Marine Information Acquisition and Security (Harbin Engineering University), Ministry of Industry and Information Technology, Harbin 150001, China
College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin 150001, China

 

† Corresponding author. E-mail: liujuan@hrbeu.edu.cn

Abstract

The multimodal admittance method and its improvement are presented to deal with various aspects in underwater acoustics, mostly for the sound propagation in inhomogeneous waveguides with sound-speed profiles, arbitrary-shaped liquid-like scatterers, and range-dependent environments. In all cases, the propagation problem governed by the Helmholtz equation is transformed into initial value problems of two coupled first-order evolution equations with respect to the modal components of field quantities (sound pressure and its derivative), by projecting the Helmholtz equation on a constructed orthogonal and complete local basis. The admittance matrix, which is the modal representation of Direchlet-to-Neumann operator, is introduced to compute the first-order evolution equations with no numerical instability caused by evanescent modes. The fourth-order Magnus scheme is used for the numerical integration of differential equations in the numerical implementation. The numerical experiments of sound field in underwater inhomogeneous waveguides generated by point sources are performed. Besides, the numerical results computed by simulation software COMSOL Multiphysics are given to validate the correction of the multimodal admittance method. It is shown that the multimodal admittance method is an efficient and stable numerical method to solve the wave propagation problem in inhomogeneous underwater waveguides with sound-speed profiles, liquid-like scatterers, and range-dependent environments. The extension of the method to more complicated waveguides such as horizontally stratified waveguides is available.

1. Introduction

Sound propagation in shallow water waveguides is an important subject in underwater acoustics. It has attracted much attention, both theoretically and experimentally, for several decades. However, the quantitative understanding of sound propagation in underwater waveguides is not available in a satisfactory way, because of the complexity of the underwater environments. For instance, sound-speed profiles as generic characteristics, and other inhomogeneities caused by organisms, vehicles, seamounts, and surface fluctuations, etc. Ray theory is one of the most utilized techniques, whose objective is tracing ray paths and thereby analyzing wave propagation, radiation, and scattering properties.[13] However, its application is limited to high frequencies, and the ray trajectories are quite sensitive to details of the media and the numerical discretization, particularly when one traces rays to long range. Normal mode method is often used in the study of underwater sound propagation in range-independent waveguides.[1,4,5] When the waveguide is range-dependent, applications of mode method is not evident. Some other classical techniques such as the finite element methods[6,7] and the finite difference methods[8,9] can be alternative to calculate wave propagation in waveguides with sound-speed profiles and inhomogeneities. Nevertheless, such methods are inconvenient when the typical wavelength is much less than the height of the waveguide or the propagation range is much larger than the height, because these circumstances require a sufficient number of calculating nodes which leads to a prohibitively expensive cost. An efficient and stable numerical method to compute the sound field in inhomogeneous underwater waveguides is significant to be addressed.

In this paper, we present the multimodal admittance method to treat the sound propagation problem in inhomogeneous underwater waveguides. The multimodal admittance method was first proposed by Pagneux et al. to compute wave propagation in a varying cross-section waveguide.[10,11] Then it has been used to analyze Lamb wave propagation in inhomogeneous elastic waveguides,[12] wave propagation in two-dimensional (2D) and three-dimensional (3D) rigid bends,[13,14] and in waveguides with penetrable scatterers.[15] In order to improve the calculation convergence, the improved version of multimodal method was developed to analyze wave propagation in waveguides with impedance boundary condition[16] or with varying cross section.[17] The multimodal admittance method stems from the multimodal method, the basic concept of which is the generalized Fourier series expansion, in other words, the solutions of the Helmholtz equation can be expanded by an orthogonal and complete basis consisting of transverse local modes. The modes are supposed to include all propagating modes and several evanescent mode. By projecting the Helmholtz equation onto the mode basis, two couple first-order differential equations with respect to the modal components of the sound pressure and its derivative are obtained. A sound condition and a radiation condition are crucial to give the initial conditions to solve the coupled mode equations. However, directly integrating the first-order differential equation from the radiation condition can bring numerically instability, because the evanescent waves are exponentially increasing in that way. The admittance matrix is introduced to evade the numerical divergence difficulty. The admittance matrix is often referred to the modal representation of the Direchlet-to-Neumann operator, linking the sound pressure and its derivative in modal representation. It can efficiently avoid the contamination generated by exponentially diverging evanescent modes. Although the multimodal admittance method has been utilized in various waveguides, the wave propagation and scattering in underwater waveguides under the effect of a sound-speed profile, still need to be investigated. For the sake of comprehensiveness, the boundary conditions of the physical problem in this paper are formulated as acoustically soft at the top water–air interface and rigid at the bottom water–rock interface, respectively. The generalization to more complicated cases, such as the ocean geometry is modeled as a two-layered waveguide, is also possible. Besides the sound-speed profiles, two inhomogeneous cases in underwater waveguides are considered: (i) waveguides with liquid-like scatters and (ii) waveguides with range-dependent environments. For the case (ii), the expansion of sound pressure by local normal modes performs poor convergence with respect to the truncation number of the transverse modes, because the local normal modes do not satisfy the natural boundary conditions of the varying cross section waveguides. The improved multimodal admittance method is presented to give a better convergence. The idea is choosing a boundary mode[18] to enrich the mode basis such that the natural boundary condition is automatically satisfied. We will show that the advantage of the improved method is that the calculation follows the classic procedure.

The paper is organized as follows: in Section 2, a brief revisiting of multimodal method is present for the wave propagation in a 2D axial-uniform waveguide with a sound-speed profile. In Section 3, we treat the wave propagation problem in a 2D waveguide with liquid-like scatterers and a soundspeed profile. Basic equations of the multimodal admittance method are given. To numerically solve the resulting coupled first-order wave equation and Riccati equation, we employ a fourth-order Magnus numerical integration scheme, which is suitable for computing wave propagation in large ranges or high frequencies. For the validation of the multimodal admittance method, numerical results of the sound fields are shown in Subsection 3.3, including a comparison of the results calculated by simulation software COMSOL Multiphysics to the results using the multimodal admittance method. In Section 4, the improved multimodal admittance method is applied to waveguides with range-dependent environments in order to realize a better convergence. The corresponding sound field is illustrated as well. Finally, in Section 5 we give the conclusions.

2. Sound propagation in waveguides with soundspeed profiles using multimodal method

The sound propagation in a 2D water-filled waveguide of semi-infinite length and finite depth h with a sound-speed profile is considered as shown in Fig. 1. The sound speed c(y) is assumed to depend merely on the depth. The boundaries are acoustically soft at y = 0 and acoustically rigid at y = h, respectively. A source is emitted from left x = 0.

Fig. 1. Configuration of a waveguide having a sound speed profile.

The harmonic time dependence in this paper is e−iωt and will be omitted in the following. The acoustic pressure p satisfies the Helmholtz equation

where ∇2 = 2/∂x2 + 2/∂y2, together with the boundary conditions

The waveguide sustains an infinite number of modes, which are complicated to be analytically computed due to the presence of sound speed gradient. We use the multimodal method to treat this problem. The basic concept of the multimodal method, based on the generalized Fourier series expansion, is to expand the unknown by a set of orthogonal and complete transverse functions (a mode basis) and the expansion coefficients (modal components) are evaluated by projecting the Helmholtz equation onto the basis.

The basis we choose obeys the eigenproblem

Here Ψn(y) denote the normal modes in a waveguide with the same dimension and boundary conditions as the model shown in Fig. 1, but having a constant sound speed . It follows

The orthogonality relation is given as

The sound pressure at any position can be expanded by Ψn(y), i.e.

where p and Ψ are vectors with elements pn and Ψn, respectively. N is the truncation number. Then projecting Eq. (1) on the basis Ψ* and taking advantage of the orthogonality relation, we get a second-order differential equation

One can also write the above equation as coupled mode equations in terms of s = xp

where , I is an identity matrix, and K2 is the ‘wavenumber’ matrix

here is a diagonal matrix with on the diagonal. The eigenvalues of K2 are exactly identical to where kx are the horizontal wavenumbers of the modes.

The matrix H can be regarded as the evolution propagator of the axial-uniform waveguide with sound-speed profiles. H is independent of x which enables us to directly compute the modal components p

where p0 is the modal components of the source condition p(x = 0), and K is the square root of K2. The choice of the sign of K is determined by the radiation condition stating that the energy radiates outwards as propagation. Since the time dependence is e−iωt, the eigenvalues of K are chosen having zero or positive imaginary part, corresponding to propagative or evanescent waves, respectively. The natural modes in the waveguide can be obtained by mutiplying the transpose of the eigenvectors of K and the modal basis Ψ. In the numerically implementation we employ a numerical integration tool called Clenshaw–Curtis interpolatory quadrature[18] to compute the integral in Eq. (8). The details of the Clenshaw–Curtis interpolatory quadrature are given in Appendix A.

Note that to solve the coupled first-order differential equations Eq. (7), originally, two initial conditions are necessary. However, due to the axial uniformity, the radiation condition is implicitly contained in the coupled equations, expressed as s/p = iK. Consequently, once the source condition p(x = 0) (the pressure distribution at x = 0) is given, the sound pressure field can be acquired by computing the modal components p (Eq. (9) and then substituting it into Eq. (5).

3. Sound propagation in waveguides with soundspeed profiles and liquid-like scatterers

In this section we investigate the sound propagation in inhomogeneous waveguides with scatterers of liquid-like material and still, depth-dependent sound speed profiles. The resultant field contains a backscattered field and a transmitted field. Using the multimodal method as mentioned in Section 2, the Helmholtz equation can be reduced to two coupled first-order differential equations. Two initial conditions, i.e., source condition and radiation condition, are crucial to solve the equations. However, numerical instability emerges when integrating the first-order differential equation from the radiation equation since the evanescent modes are exponential growing in that way (from right to left). The problem can be circumvented by applying the admittance matrix.

3.1. Multimodal admittance method

The configuration is illustrated in Fig. 2. Ωs and Ωres denote the space occupied by the scatterers and the exterior space outside the scatterers in the waveguide, respectively. All the scatterers are involved in the region x ∈ [0,xmax]. Therefore the resultant scattered field contains a backscattered field in the region 0 ⩽ xxmax and a transmitted field in the region xxmax. The number and the shape of the scatterers are arbitrary. Without loss of generality, we consider only one scatterer in the following. The sound speed and mass density of the scatterer are c0(y) and ρ0, respectively, which are in contrast to the sound speed c(y) and density ρ for the host medium. ρ0 and ρ are constant. The sound pressure p(x,y) in such an inhomogeneous waveguide satisfies the Helmholtz equation

Fig. 2. Configuration of an inhomogeneous waveguide with a sound-speed profile and a liquid-like scatterer.

where 1 + α(y) = c(y)/c0(y), and in case of absorbing scatterers, α(y) is a complex function with ℑm(α) > 0. The boundary conditions are imposed to be acoustically rigid for the bottom wall and soft for the top wall, i.e.,

and the continuous conditions on the interface of the scatterer and the host medium are given as

where n represents the exterior normal of the scatterers and p+ (p) denotes the pressure at the exterior (interior) side of the interface of the scatterer and the host medium.

First we rewrite Eq. (10) into weak formulation form by introducing a test function q(x,y)[15] and integrating it in the whole region. This leads to

where q(x,y) = Q(x)Ψn(y), Q(x) is compactly supported, and ξ = (ρ/ρ0)(1+α)2 − 1. Ω = ΩsΩres is the whole region of the waveguide. We expand p by the basis Ψn(y) (Eq. (3))

where p is a vector with pn as the elements and N is the truncation number. Then substituting Eq. (14) into Eq. (13), after some manipulations, we obtain

where

Here a(x) and b(x) are the shape functions of the scatterer (see Fig. 2). The integrals in E and the left term of D can be analytically computed a(x)b(x)ψnψmdy=yh(sinc((mn)πyh)sinc((m+n+1)πyh))|y=a(x)y=b(x),a(x)b(x)ψnyψmydy=π2(m+12)(n+12)yh3×(sinc((mn)πyh)+sinc((m+n+1)πyh))|y=a(x)y=b(x).

K2 and the rest term of D are numerically calculated by the Clenshaw–Curtis interpolatory quadrature given in Appendix A. When there is no scatterer, i.e. a(x) = b(x), ρ = ρ0, and α = 0, such that E = I and D = 0, the problem turns to be the wave propagation problem described by Eq. (7). The second-order differential equation (Eq. (15)) is actually the reformulation of the Helmholtz equation Eq. (10) via projecting it onto the constructed basis Ψn(y).

We write Eq. (15) into two coupled first-order evolution equations by assuming s = E∂xp

where , and H can be considered as the evolution propagator of the corresponding waveguide. The propagator here is x-dependent while the propagator of the axial-uniform waveguide (Eq. (7)) is x-independent. Two initial conditions are required to solve the equations, a source condition at x = 0 and a radiation condition at xxmax. Note that equation (16) cannot be integrated directly because the evanescent modes are exponentially increasing when integrating from right to left which leads to numerical divergence.[12,19] Nevertheless, evanescent modes are not negligible since they are significant for the near field of any inhomogeneity in the waveguide. We employ the admittance matrix to avoid this circumstance. The admittance matrix Y (x) is the modal representation of Direchlet-to-Neumann operator which is defined as s = Y p, linking the acoustic pressure p to its x-derivative s.

Substituting s = Y p into Eq. (16), we yield two equations: the Riccati equation with respect to Y and a first-order differential equation governing the modal components p

Given an initial value of Y, the admittance matrix for all x is able to be obtained. The initial value of Y can be extracted from the transmitted field region xxmax where there only exist right-going waves with E(xxmax) = I and D(xxmax)=0. The initial condition follows

where K is the square root of K2. The values of K are chosen to have eigenvalues with zero or positive imaginary part. Therefore, it is available to compute p once the source condition is given. Note that the required source condition is the sound pressure distribution at x = 0, meaning that it is the superposition of the incident wave and the reflected wave. The modal components of the source condition can be obtained with the aid of the reflection matrix R

where R = (iK + Y )−1(iKY) and pinc is the modal components of the incident wave.

3.2. Numerical integration

The differential equations (17a) and (17b) cannot be readily solved because the matrices E, D, and Y depend on x. We give the numerical procedure to solve these equations, and one can follow the same steps to compute the wave field governed by Eqs. (40) and (41) in Section 4. We assume that an operator G is defined as p(xmax) = G(xmax,x)p(x) for xxmax. Thus G satisfies the first-order differential equation xG = −GE−1Y. G can be numerically computed by integrating from right at x = xmax to left at x = 0 with the initial condition G(xmax, xmax) = I. The modal components p can be obtained thereafter.

Inspired by Refs. [20,21], we use the fourth-order Magnus method for numerical integration of linear differential equations. The fourth-order Magnus method is suitable for the situation where the calculating range is much larger than the typical wavelength, or the calculated propagation distance is much larger than the height of the waveguide, since it necessitates less calculating nodes than several classical numerical techniques, such as finite difference methods and finite element methods.[20,22] The detail of this numerical scheme to calculate Y (x), G(x), and p is now presented. The interval x ∈ [0, xmax] is discretized by a series of coordinates x1 > x2 > ··· > xN with xN = 0 and x1 = xmax. Equation (16) can be rewritten as

Here the Magnus matrix Mn is given as

where

are the evaluations of the matrix H at the nodes of the fourthorder Gauss–Lagendre quadrature in the segment from xn to xn+1, and δn = xn+1xn. This iteration is operated from right to left such that δn < 0. Note that when Mn is x-independent, which implies that there is no scatterer in the medium, the results derived from the fourth-order Magnus scheme are exact. The matrix exponential eMn can be devided as

Substituting Eq. (21) into Eq. (20) and taking using of s = Y p, Y (x) and G(x) then can be derived as

Therefore p is able to be solved benefiting from the operator G(x) since G(x) establishes the relation between p(xi) and p(xj) at any two positions in the range [0, xmax]. The iteration result for p is given as

The iteration of p is operated from 0 to xmax, starting from the modal components p0 of source condition.

To summarize, once the geometries and characteristic parameters (sound speed and density) of the host medium and liquid-like scatterers are known, Y (x) and G(x) in the range [0, xmax] are determined by integrating from the radiation condition. And the sound field generated by any possible source in this waveguide can be obtained by iteration from source condition. In numerical implementation, the series of Eq. (14) is truncated at finite N. The choice of N needs to make sure all the propagative modes and several evanescent modes are included in the series.

3.3. Numerical results and method validation

Normally, a point source is appropriate to use in practical problems in underwater acoustics. Assuming a point source is located at (x,y) = (0, ys), we give the expanding expression of the corresponding Green’s function G(x,y) on the constructed basis Ψn(y) (more details see Appendix B)

where g is the modal component vector g ≡ (gn)

Figure 3 displays the examples of sound propagation in inhomogeneous waveguides with or without scatterers, generated by a point source at y = 0.2h, where h = 30 m is the depth of the waveguides and the frequency is f = 500 Hz. Figures 3(a) and 3(b) show the sound field of a waveguide having negative sound-speed profile gradient c = 1500 × (1 − 10−3y) and in the absence of scatterer. The calculation in Fig. 3(a) has been done using the multimodal method as mentioned in Section 2, while using a commercial software COMSOL Multiphysics in Fig. 3(b). It is observed that the results agree well. The acoustic wave bends locally towards the bottom boundary where the sound speed is low, and then reflected by the surfaces. The waveguide considered in Figs. 3(c) and 3(d) contains a scatterer formulated by , and parameterized by ρ/ρ0 = 0.3 and α = 0.5, and the sound speed profile in the host medium is still c = 1500 × (1 − 10−3y). The sound fields illustrated in Figs. 3(c) and 3(d) are computed by the multimodal admittance method and commercial software COMSOL Multiphysics, respectively. We can observe that the result in Fig. 3(c) is extremely close to that in Fig. 3(d). Moreover, the change of the wavelength between the host medium and the scatterer is clear. It is verified that the multimodal admittance method is a powerful and efficient way to analyze wave propagation in inhomogeneous waveguides with sound-speed profiles.

Fig. 3. Sound fields (pressure amplitude) in inhomogeneous waveguides having sound-speed profiles without ((a) and (b)) and with ((c) and (d)) a scatterer. The sound source is a point source located at y = 0.2h where h = 30 m is the depth. The frequency is f = 500 Hz. The sound-speed profiles in the host medium are assumed to be the same in both cases expressed as c = 1500 × (1 − 10−3y) as shown in the left panel. (a) Sound field in the waveguide with no scatterer. The result is computed by the multimodal method; (b) Sound field in the waveguide with no scatterer. The result is computed by a numerical calculating software COMSOL Multiphysics; (c) Sound field in the waveguide containing a scatterer with parameters ρ/ρ0 = 0.3 and α = 0.5. The result is calculated by multimodal admittance method; (d) Sound field in the waveguide with the same scatterer as in panel (c). The result is calculated by COMSOL Multiphysics.

Figure 4 depicts the scattering field of a waveguide with an absorbent scatterer, excited by a point source located at y = 0.2h, where h = 30 m is the depth. The frequency is f = 500 Hz. The sound speed of the host medium is same as that in Fig. 3, expressed as c = 1500 × (1 − 10−3y). The parameters of the scatterer are set as ρ/ρ0 = 0.3 and α = 0.01i. The complex sound speed in the scatterer means that it is lossy. The results are computed by the multimodal admittance method. When sound propagates across the absorbent scatterer, sound interacts with the scatterer and the wavelength changes in the scatterer range. Figure 4(b) exhibits the corresponding variation of energy flux as a function of x of Fig. 4(a). In inhomogeneous waveguides the expression of the energy flux is written as

Fig. 4. Sound propagation in a waveguide with a sound-speed profile and an absorbent scatterer: ρ/ρ0 = 0.3 and α = 0.01 i, generated by a point source. The source is located at 0.2h, where h = 30 m. The frequency is f = 500 Hz. The sound speed of the host medium is c = 1500 × (1 − 10−3y). (a) Sound field (pressure amplitude); (b) Energy flux as a function of x of the waveguide in panel (a).

and the reciprocity an energy conservation have been proved in Refs. [15, 23]. It is shown that the sound energy suffers loss when propagation across the absorbent scatterer, and sound propagates without loss in other regions due to the perfectly total reflection boundaries.

4. Sound propagation in waveguides with range-dependent environments and sound-speed profiles
4.1. The improved multimodal admittance method

In the shallow ocean, the top and bottom surfaces are important since sound inevitably bounces by the surfaces. There are some situations where the cross section of the waveguide is range-dependent, and therefore strongly influence the sound field. In the section, the wave propagation in a waveguide with continuously and varying cross section is inspected. For simplicity, we only formulate the waveguide with varying bottom wall in the region x0xxmax (see Fig. 5). The top wall is acoustically soft, and the bottom is rigid, described by h(x). A stable source is given at x = 0. The sound speed depends on the depth. The corresponding sound pressure satisfies the Helmholtz equation

Fig. 5. Configuration of a varying cross-section waveguide with sound-speed profiles.

with the boundary conditions

where

First we construct a mode basis ψn(y;x) to express p and s. Here ψn(y;x) are local normal modes, which correspond to the eigenfunctions of the eigenproblem with homogeneous Direchlete boundary condition for top wall ψn(0;x) = 0 and Neumann boundary condition for bottom wall ∂ψn(h(x);x)/∂y = 0 at each local x. It follows

Note that the natural boundary condition imposed on the bottom surface is ∂p(h(x);x)/∂y = h′(x)∂p(h(x);x)/∂x, but the basis functions ψn(y;x) do not satisfy it such that the series representation of p and s onto the basis ψn(y;x) (n = 0,1,2,…,N−1) will exhibit poor convergence properties. We use the improved version of the multimodal admittance method[17] to treat this problem. The idea is adding a supplementary function ψ−1, so-called a boundary mode, to enrich the mode basis. It is easy to write such a boundary mode ψ−1(y;x) by Schmidt orthogonalization

where χ satisfies ∂χ(h(x);x)/∂y ≠ 0 and χ(0;x) = 0. χ is arbitrary, and here we choose

where a−1 is the normalization factor to make sure ψ−1 normalized. Therefore we get a basis ψn(n = −1,0,1,…,N−1) with orthogonality relation

Then we write the Helmholtz equation into the first-order evolution form in terms of s = xp

p and s can be written as the expansion on ψn(y; x)

The expansion convergence increases from 1/N2 to 1/N4 by adding a boundary mode (see Appendix C). Substituting Eq. (35) into Eq. (34) and integrating over the depth of the waveguide, we yield

Write Eq. (36) into two coupled first-order differential equations governing the modal components p ≡ (pn) and s ≡ (sn)

where the evolution propagator , and the elements in A and K2 are

The explicit expressions for the formulas are given in Appendix D. Note that for 0⩽ n, m ⩽ (N − 1),

correspond to the square of the usual horizontal wavenumbers of the modes. For n = m = −1,

is found negative, corresponding to an evanescent mode, such that the boundary mode ψ−1 automatically satisfies the radiation condition (for details see Ref. [17]).

Equation (39) implies that the Helmholtz equation is posed as initial value problems of the coupled first-order differential equations. Equation (38) cannot be integrated directly, similar to the discussion in Section 3. We introduce the admittance matrix Y as

Inserting Eq. (39) into Eq. (38), it follows a Riccati equation with respect to Y

One can get the admittance matrix Y at any x with an initial condition. The initial condition of Y can be extracted from the radiation condition in the constant cross-section region xxmax

p eventually satisfies the first-order evolution equation

The sound field can be computed once the source condition p(x = 0) is given. We note that the properties of the boundary mode ψ−1, orthogonal to the rest N modes and evanescent as propagation, lead to the classic form of coupled first-order evolution equations. The numerical regime to solve Y and p follows the same steps as described in Subsection 3.2.

To summarize, using the (improved) multimodal admittance method, the wave propagation problem for an inhomogeneous waveguide is effectively reduced to initial value problems of two coupled first-order evolution equations with respect to the modal components of sound field. The evolution behavior is well described by an evolution propagator H. The inhomogeneities influence strongly on the evolution propagator. For the axially uniform waveguide, the evolution propagator (Eq. (7)) is apparently x-independent. The effect of scatterers lies on the anti-diagonal block (Eq. (16)), and the range-dependent cross-section leads to additive terms on the main diagonal block (Eq. (38)). It makes the extension of the multimodal admittance method to the wave propagation in the underwater waveguides with both liquid-like scatterers and varying cross-section available in a more straightforward way.

4.2. Numerical results

The sound field in a varying cross section waveguide, generated by a point source located at 0.2h, h = 50m, is plotted in Fig. 6. The frequency is f = 500 Hz. The waveguide has a negative sound-speed profile c = 1500 × (1−10−3y). There is a seamount at the bottom surface, described as h(x) = hh/10(1 − cos(2π (xh)/h)), h < x < 2h, and h = 50 m elsewhere. The truncation number is N = 200, in order to guarantee that all the propagative modes and a few evanescent modes are included. It is apparent that the sound interact with the seamount and suffer strong reflection. Figure 6 shows the energy flux ratio at any position x with respect to the energy flux at x = 0, i.e., Ef(x)/Ef(0). The energy is conserved.

Fig. 6. (a) Sound propagation in a waveguide with a sound-speed profile and a seamount, generated by a point source. The source is located at 0.2h, where h = 50 m. The frequency is f = 500 Hz. The sound speed is c = 1500 × (1 − 10−3y). (b) Energy flux ratio Ef(x)/Ef(0) as a function of x for the waveguide in panel (a).
5. Conclusion

We have presented the multimodal admittance method to analyze the sound propagation in underwater waveguides. Various aspects of inhomogeneity in ocean acoustics have been taken into account: sound speed profiles, liquid-like scatterers and range-dependent bathymetry. The improved multimodal admittance method by adding a boundary mode on the mode basis is inspected to realize better convergence in the varying cross section cases. We show that using the multimodal admittance method, the Helmholtz equation is transformed into two first-order coupled differential equations with respect to the modal components of sound pressure and its x-derivative on the constructed mode basis. These first-order coupled differential equations govern the sound evolution behavior for the specified inhomogeneous waveguide through an evolution propagator. Diverse inhomogeneity has different effect on the evolution propagator. For the axially uniform waveguide with only depth-dependent sound speed profiles, the evolution propagator is x-independent. When penetrable scatterers are imposed in the waveguide, the evolution propagator will add extra terms on the anti-diagonal block. The range-dependent cross-section will lead additive terms on the main diagonal block.

The series expansion with respect to transverse modes contains evanescent modes. The evanescent modes can cause numerical instability during the integration from radiation condition. The admittance matrix is introduced to stably compute the integrals, and it is able to eliminate the contamination of exponentially diverging of evanescent modes. The admittance matrix links the modal components of sound pressure and its x-derivative on the constructed mode basis. To summary, the key point of the multimodal admittance method is projecting the Helmholtz equation onto the mode basis and using the admittance matrix to compute the solutions. The multimodal admittance method has significant advantages compared with some typical methods applied in underwater acoustics. There is no assumption about frequencies in the derivation of the multimodal admittance method, which enables us compute the sound field excited by sources with arbitrary frequency. Besides, the multimodal admittance method necessitates less calculating nodes than several classical numerical techniques, such as finite difference methods and finite element methods, which reduce the computational expense.

This method is able to be applied to more complicated cases where the modes in the waveguides are not orthogonal. We can construct a basis and make the projection based on the bi-orthogonality of the local mode basis. The method gives the possibility to analyze the physical systems no matter they are hermitian or not. In the paper, we choose 2D planar waveguides with simple Dirichlet and Neumann boundary conditions as models to demonstrate the multimodal admittance method. The simple models are easy to be manipulated and accessed. All the formulations presented in the paper are analytical, the existing errors result from the use of truncation number (the projection) and the numerical integrals. The extension to more complicated situations, such as three-dimensional (3D) cylindrical waveguides and horizontally layer waveguides, is also available.

Reference
[1] Jensen F B Kuperman W A Porter M B Schmidt H 2011 Computation ocean acoustics 2 New York American Institute of Physics 12 15 https://link.springer.com/book/10.1007/978-1-4419-8678-8
[2] ČErvený V 2001 Seismic Ray Theory Cambridge University Press 1098 1108 https://doi.org/10.1007/0-387-30752-4_134
[3] Štumpf M De Hoop A T Vandenbosch G A E 2013 IEEE Trans. Anten. Propag. 61 2676
[4] Yang S E 2009 Theory of underwater sound propagation Harbin Harbin Engineering University Press 6 in Chinese
[5] Mo Y X Piao S C Zhang H G Li L 2014 Acta Phys. Sin. 63 214302 in Chinese
[6] Ihlenburg F Babuska I 1997 SIAM J. Numer. Anal. 34 315
[7] Ihlenburg F Babuska I 2000 SIAM Rec. 42 451
[8] Singer I Turkel E 1998 Comput. Method Appl. M 163 343
[9] Nabavi M Siddiqui M H K Dargahi J 2007 J. Sound Vib. 307 972
[10] Pagneux V Amir A Kergomard J 1996 J. Acoust. Soc. Am. 101 2034
[11] Amir A Pagneux V Kergomard J 1997 J. Acoust. Soc. Am. 101 2504
[12] Maurel A Pagneux V 2002 P. R. Soc. A-Math. Phys. 458 1913
[13] Félix S Pagneux V 2001 J. Acoust. Soc. Am. 110 1329
[14] Félix S Pagneux V 2002 Wave Motion 36 157
[15] Maurel A Mercier J F Félix S 2014 J. Acoust. Soc. Am. 135 165
[16] Bi W P Pagneux V Lafarge D Aurégan Y 2007 J. Acoust. Soc. Am. 122 280
[17] Maurel A Mercier J F Pagneux V 2014 P. R. Soc. A-Math. Phys. 470 20130448
[18] Waldvogel J 2006 BIT 46 195
[19] Pagneux V Maurel A 2006 P. Roy. Soc. A-Math. Phys. 462 1315
[20] Lu Y Y 2005 J. Comput. Appl. Math. 172 247
[21] Iserles A Marthinsen A N?rsett S P 1999 BIT 39 281
[22] Lu Y Y McLaughlin J R 1996 J. Acoust. Soc. Am. 100 1432
[23] Pagneux V Maurel A 2004 J. Acoust. Soc. Am. 116 1913